library(gstat)
library(lattice)

d <- read.csv("C:\\temp\\gw.csv")
names(d) <- c('loc','well','y','x','z','date')
coordinates(d) <- ~ x+y
d$x <- coordinates(d)[,1]
d$y <- coordinates(d)[,2]

bbox <- summary(d)$bbox
cellsize <- 200
buffer = 1000
xmin <- bbox[1] - buffer
ymin <- bbox[2] - buffer
xrange <- bbox[3]-bbox[1]
yrange <- bbox[4]-bbox[2]
xdim <- ceiling((xrange+buffer*2)/cellsize)
ydim <- ceiling((yrange+buffer*2)/cellsize)
gt <- GridTopology(c(xmin,ymin),c(cellsize,cellsize),c(xdim,ydim))
gd <- SpatialGrid(gt)

# check for normality
densityplot(d$z)
densityplot(log(d$z))

d$logz <- log(d$z)

v <- variogram(logz ~ 1, data=d, width=300)
plot(v, plot.numbers=T)

vm1 <- vgm(psill=0.4, model="Gau",nugget=0.1, range=6000)
plot(v, model=vm1)

vm2 <- fit.variogram(v, vm1)
#plot(v, model=vm2)

g <- gstat(id="logz", formula=logz ~ 1, data=d, model=vm1)
okr <- predict(g, id="logz", newdata=gd)
okr$z.pred <- exp(okr$logz.pred)

#pdf("C:\\temp\\gw.pdf", width=17,height=11)

#plot.new()
#trellis.par.set(sp.theme())
#spplot(okr, "z.pred")
#contour(okr,"z.pred",add=T)
#plot(d, add=T, cex=0.6, pch=16)

#spplot(okr, "logz.var")

invdist <- idw(logz ~ 1, d, gd, maxdist=2000)
invdist$z.pred <- exp(invdist$var1.pred)
#spplot(invdist, "z.pred")
#contour(invdist,"z.pred",add=T)
#plot(d, add=T, cex=0.6, pch=16)

# convert to matrix
#gm <- matrix(okr$z.pred, ncol=xdim, byrow=T)
#persp(gm)

# Use RGL to display in 3D
library(rgl)
#okr <- invdist
rgl.points(coordinates(d)[,1],d$z*-80, coordinates(d)[,2], size=3,col="blue")
rgl.surface( ((1:xdim)*cellsize)+bbox[1]-buffer,
             ((ydim:1)*cellsize)+bbox[2]-buffer,
              okr$z.pred*-80)
axes3d()


